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Abstract: Over the last fifty years, the finite element method has emerged as a suc- 
£SJ ' cessful methodology for solving a wide range of partial differential equations. At the 

heart of any finite element simulation is the assembly of matrices and vectors from 
finite element variational forms. In this paper, we present a general and unified frame- 
work for finite element assembly. Based on this framework, we propose a specific 
k>( ■ software interface called Unified Form-assembly Code (UFC) between problem-specific 

' and general-purpose components of finite element programs. The interface is general 

| in the sense that it applies to a wide range of finite element problems (including mixed 

finite elements and discontinuous Galerkin methods) and may be used with libraries 
that differ widely in their design. The interface consists of a minimal set of abstract 
C++ classes and data transfer is via plain C arrays. 

We discuss how one may use the UFC interface to build a plug-and-play system for 
finite element simulation where basic components such as computational meshes, linear 
algebra and, in particular, variational form evaluation may come from different libraries 
and be used interchangeably. We further discuss how the UFC interface is used to glue 
together components from the FEniCS suite of software to provide an integrated high- 
level environment where variational forms may be entered as expressions directly in 
Python and assembled efficiently into sparse matrices. 

A central design goal for the interface is to minimize dependency on external libraries 
for the problem-specific code used in applications. Thus, the UFC interface consists 
of a single C++ header file and does not rely on external libraries for its operation. 
In particular, the UFC interface does not depend on any other FEniCS components. 
As a result, finite element code developers may use the interface to detach equation- 
specific details from general-purpose library code, allowing very flexible connections to 
alternative libraries. We encourage developers of finite element libraries to incorporate 
the interface in their libraries. The UFC interface is released into the public domain. 

Keywords: Finite elements, assembly, implementation, code generation, UFC 

Reference to this paper should be made as follows: M. S. Alnaes, A. Logg, K.-A. 
Mardal, O. Skavhaug, and H. P. Langtangen (2008) 'Unified Framework for Finite 
Element Assembly'. 



1 Introduction 



Software for solving physical problems have traditionally 
been tailored to the problem at hand, often resulting in 
computationally very efficient special-purpose codes. How- 
ever, experience has shown that such codes may be difficult 
and costly to extend to new problems. To decrease turn- 
over time from problem definition to its numerical solution, 
scientific code writers have to an increasingly larger extent 
tried to create general libraries, containing common nu- 
merical algorithms applicable to a wide range of problems. 
Such libraries can reduce the size of the application code 
dramatically and hide implementation details. In the field 
of finite element solution of partial differential equations, 
many general and successful libraries have emerged during 
the last couple of decades, e.g., Cactus, Cogito, COMSOL 
Multiphysics, Deal.II, Diffpack, DOLFIN (FEniCS), Get- 
fcm++, Kaskade, Sundance, and UG (see the reference list 
for papers and websites). 

General finite element libraries implement many stan- 
dard mathematical and numerical concepts, but the soft- 
ware components are often not as carefully designed as 
their mathematical counterparts. From a software engi- 
neering point of view it is important to achieve clear sepa- 
ration of the various software components that build up a 
finite element library, such that each component can be re- 
placed separately. Not only does this offer greater flexibil- 
ity for application and library developers, but it also makes 
the software easier to maintain, especially under changing 
requirements of several developers in long-term projects. 
These arguments have received much attention by dcvcl- 
oper s of general finit e elem e nt libraries in r e cent ye ars (see, 
e.g.. lBangerth et al.l (|2007t ): iBastian et all l|2007al )). 

Well designed libraries provide clear interfaces to repre- 
sent this separation. Typically, the application code uses 
functions or objects in the interface to perform basic "high 
level" steps of the solution process. Problem-specific de- 
tails, such as the variational form, the mesh and coefficients 
are passed through the interface to the library to compute 
a solution. Such libraries and their interfaces are generally 
referred to as problem solving environments (PSEs). 

However, one fundamental issue in designing such soft- 
ware libraries is how to separate problem-specific code from 
general library code. Some components, such as computa- 
tional meshes and linear algebra, may be implemented as 
reusable components (e.g. as a set of C++ classes) with 
well-defined interfaces. However, other components, such 
as variational forms, are intrinsically problem-specific. As 
a result, those components must either be implemented 
and provided by the user or generated automatically by 
the library from a high-level description of the variational 
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form. In cither case, it becomes important to settle on a 
well-defined interface for how the library should communi- 
cate with those problem-specific components. 

The design of such an interface is the subject of the 
present paper. We propose a C++ interface called UFC, 
which provides an interface between general reusable fi- 
nite element libraries and problem-specific code. In other 
libraries, the implementation of finite elements and vari- 
ational forms is usually tied to the specific mesh, matrix 
and vector format in use, while in UFC we have strived to 
decouple these concepts. Furthermore, the interface is de- 
signed to allow for a variety of elements such as continuous 
and discontinuous Lagrange, Nedelec and Raviart-Thomas 
elements. 

To make a successful interface, one needs a sufficiently 
general framework for the underlying mathematical struc- 
tures and operations. The software interface in the cur- 
rent paper relies on a more general view of variational 
forms and finite element assembly than commonly found in 
textbooks. We therefore precisely state the mathematical 
background and notation in Sections [5] and [3J 

UFC is significantly inspired by our needs in the 
tools FFC, SFC, and DOLFIN, which are softwar e 
units within FEniCS, se e iKirbv and LogsJ (|2006l 120071 ); 
Alnaes and Mardall (|2008l ) ; lLogsJ (|2007l) . The interplay be- 
tween these tools and UFC is explained in Section [H which 
provides additional and more specific motivation for the 
design of UFC. Highlights of the interface are covered in 
Section \5[ Section [5] contains some examples of high-level 
specifications of variational forms with the form compilers 
FFC and SFC, which automatically generate code compat- 
ible with the UFC interface for computing element matri- 
ces and vectors. We also explain how the interface can be 
used with existing libraries. 

We here note that as a result of the UFC interface, the 
two form compilers FFC and SFC may now be used in- 
terchangeably since both generate code conforming to the 
UFC interface. These form compilers were developed sep- 
arately and independently. The work on UFC was initially 
inspired by our efforts to unify the interfaces for these form 
compilers. 

Related Work 

One major reason for the success of general finite element 
libraries is that many widely different physical problems 
can be solved by quite short application codes utilizing 
the same library. The opposite strategy, i.e., one appli- 
cation utilizing different alternative libraries, has received 
less attention. For example, an application might want 
to use an adaptive mesh data structure and its function- 
ality from one library, a very efficient assembly routine 
from another library, basic iterative methods from, e.g., 
PETSc, combined with a preconditioner from Trilinos or 
Hypre. To make this composition a true plug-and-play op- 
eration, the various libraries would need to conform to a 
unified interface to the basic operations needed in finite 
element solvers. Alternatively, low level interfaces can be 



implemented with thin wrapper code to connect separate 
software components. 

In numerical linear algebra, the BLAS and LAPACK in- 
terfaces have greatly simplified code writing. By express- 
ing operations in the application code in terms of BLAS 
and LAPACK calls, and using the associated data (array) 
formats, one program can be linked to different implemen- 
tations of the BLAS and LAPACK operations. Despite 
the great success of this approach, the idea has to little ex- 
tent been explored in other area s of computationa l science. 
One recent example is Easy viz ( Ring et all ( 2008 )). a thin 
unified interface to curve plotting and 2D/3D scalar- and 
vector-field visualization. This interface allows an applica- 
tion program to use a MATLAB-compatible syntax to cre- 
ate graphics, independently of the choice of graphics pack- 
age (Gnuplot, Grace, MATLAB, VTK, Visit Open DX, 



2 Finite Element Discretization 



etc.). Another example is GLAS (|Meerbergenj ( 20081 )). a 
community initiative to specify a general interface for lin- 
ear algebra libraries. GLAS can be viewed as an extension 
and modernization of the BLAS/LAPACK idea, utilizing 

powerful constructs offered by C ++. 

Within finite elements, DUNE (jBastian et al. I (j2007allbh 
is a very promising attempt to define unified interfaces 
between application code and libraries for finite clement 
computing. DUNE provides interfaces to data structures 
and solution algorithms, especially finite element meshes 
and iterative solution methods for linear systems. In prin- 
ciple, one can write an application code independently of 
the mesh data structure and the matrix solution method. 
DUNE does not directly address interfaces between the 
finite element problem definition (element matrices and 
vectors), and the assembly process, which is the topic of 
the present paper. Another difference between DUNE and 
our UFC interface is the choice of programming technol- 
ogy used in the interface: DUNE relies heavily on inlining 
via C++ templates for efficient single-point data retrieval, 
while UFC applies pointers to chunks of data. However, 
our view of a finite element mesh can easily be adapted 
to the DUNE-Grid interface. The DUNE-FEM module 
(under development) represents interfaces to various dis- 
cretization operators and serves some of the purposes of 
the UFC interface, though being technically quite differ- 
ent. 

In the finite element world, there are many competing 
libraries, each with their own specialties. Thin interfaces 
offering only the least common denominator functional- 
ity do not support special features for special problems 
and may therefore be met with criticism. Thick interfaces, 
trying to incorporate "all" functionality in "all" libraries, 
become too complicated to serve the original purpose of 
simplifying software development. Obtaining community 
consensus for the thickness and syntax of a unified inter- 
face is obviously an extremely challenging process. The 
authors of this paper suggest another approach: a small 
group of people defines a thin (and hence efficient and easy- 
to-use) interface, they make the software publicly available 
together with a detailed documentation, and demonstrate 
its advantages. This is our aim with the present paper. 



2.1 The Finite Element 

A finite element is mathematically defined as a triplet con- 
sisting of a polygon, a polyn omial function space, and a 
set of linear functionals, see ICiarlet (Il978l) . Given that 
the dimension of the function space and the number of the 
(linearly independent) linear functionals are equal, the fi- 
nite element is uniquely defined. Hence, we will refer to a 
finite element as a collection of 

• a polygon K, 

• a polynomial space Vk on K, 

• a set of linearly independent linear functionals, the 
degrees of freedom, £i : Vk — ► R, i = 1, 2, . . . , uk ■ 

With this definition the basis functions {<pf}2fi are ob- 
tained by solving the following system of equations, 



,n K . 



(1) 



The computation of such a nodal basis can be auto- 
mated, given (a basis for) the polynomial space Vk an d 
the set of linear functio nals {^i}™^, see iKirbvl (|2004l ); 
Alnaes and Mardall (|2008h . 



2.2 Variational Forms 

Consider the Poisson problem - V-(wVu) = / with Dirich- 
let boundary conditions on a domain fi C R d . Multiplying 
by a test function v G Vh and integrating by parts, one 
obtains the variational problem 



wVw • Vwfc, dx 



vfdx, \fv G Vh 



(2) 



for the approximation Uh G Vh- If iw, / € Wh for some 
discrete function spac^H Wh we may thus write @ as 

a(v,u h ;w) = L(v; f) Vv G Vh, (3) 

where the trilinear form a : Vh, X Vh X Wh — > R is given by 

(4) 



(5) 



a(v,Uh',w) = / wVd • Vuh Ax 
Jn 

and the bilinear form L : Vh x Wh — > i? is given by 



L(v;f) 



vf dx. 



Note here that a is bilinear for any given fixed w G Wh 
and L is linear for any given fixed / G Wh- 
in general, we shall be concerned with the discretization 
of finite element variational forms of general arity r+n > 0, 

a-.VjxVZx-.-x VI x Wl x W 2 h x ■ • ■ x W h ' -> M, (6) 



x It is assumed that any given function may be represented (ex- 
actly or approximately) in some finite element space. Alternatively, 
functions may be approximated by quadrature. Quadrature represen- 
tation is not discussed here, but is covered by the UFC specification 
and implemented by the form compilers FFC and SFC. 



defined on the product space X V? X • • • X V£ X x Wfi x 
• • • x W£ of two sets {VI }j =1 , {Wl }™ =1 of discrete function 
spaces on f2. We refer to (v\, V2, ■ ■ ■ , v r ) G x V 2 x • • • x V£ 
as primary arguments, and to (w%, W2, ■ • ■ , w n ) G x 
W% x • • • x as coefficients and write 

fl = a(«i,...,w r ;wi,...,w„). (7) 

In the simplest case, all function spaces are equal but there 
are many important examples, such as mixed methods, 
where the arguments come from different function spaces. 
The choice of coefficient function spaces depends on the 
application; a polynomial basis simplifies exact integration, 
while in some cases evaluating coefficients in quadrature 
points may be required. 

2.3 Discretization 

To discretize the form a, we introduce a set of 
bases {<!>]} f =1 , {$ }£L 1; • . . , {4>\}f=i f° r the function spaces 
Vfr,Vfi, . . . , V£ respectively and let i = (ii, is, . . . , i r ) be a 
multiindex of length \i\ = r. The form a then defines a 
rank r tensor given by 

A l =a((j)l i ,$ 2 ,...,(t) r tr ;w 1 ,W2,-.-,w n ) Mi el, (8) 

where 1 is the index set 

{(1,1,...,1),(1,1,...,2),...,(A\A 2 ,...,AH}. 

We refer to the tensor A as the discrete operator generated 
by the form a and the particular choice of basis functions. 
For any given form of arity r+n, the tensor A is a (typically 
sparse) tensor of rank r and dimension IV^I x | V h 2 | x . . . x 
\V£\ = N 1 x N 2 x ... x N r . 

Typically, the rank r is 0, 1, or 2. When r = 0, the 
tensor A is a scalar (a tensor of rank zero), when r = 1, 
the tensor A is a vector (the "load vector" ) and when r = 2, 
the tensor A is a matrix (the "stiffness matrix" ) . Forms of 
higher rank also appear, though they are rarely assembled 
as a higher-dimensional sparse tensor. 

Note here that we consider the functions wi, W2, • • • , w n 
as fixed in the sense that the discrete operator A is com- 
puted for a given set of functions, which we refer to as 
coefficients. As an example, consider again the varia- 
tional problem ([2]) for Poisson's equation. For the trilinear 
form a, the rank is r = 2 and the number of coefficients is 
n = 1, while for the linear form L, the rank is r = 1 and 
the number of coefficients is n = 1. We may also choose 
to directly compute the action of the form a obtained by 
assembling a vector from the form 

a{v\;w\,w2) = / w±Vvi ■ V1U2 da;, (10) 
Jn 

where now r — 1 and n = 2. 

We list below a few other examples to illustrate the no- 
tation. 



Example 2.1. Our first example is related to the diver- 
gence constraint in fluid flow. Let the form a be given by 

a(q,u) = / qV-udx, q(=V£, ueV, 2 , (11) 

JQ. 

where is a space of scalar-valued functions and where 
V 2 is a space of vector-valued functions. The form a : 

V h x V h M has two 

primary arguments and thus r = 2. 
Furthermore, the form does not depend on any coefficients 
and thus n = 0. 

Example 2.2. Another common form in fluid flow (with 
variable density) is 

a(v,u;w,g)= / g (w ■ Vu) ■ v dx. (12) 
Jn 

Here, v € V^, u e V 2 , w € W\, g e W 2 , where V } ] , V 2 , 
and are spaces of vector-valued functions, while W 2 
is a space of scalar-valued functions. The form takes four 
arguments, where two of the arguments are coefficients, 

a-.V^x V 2 x Wl x W 2 E. (13) 

Hence, r = 2 and n = 2. 

Example 2.3. We next consider the following form ap- 
pearing in nonlinear convection- diffusion with a power-law 
viscosity, 

a(v;w, n, g) = I g{w- Vw) • v + Vw| 2<? Vu> : Vv dx. (14) 
Jn 

Here, v € V£, w € W\, fi G W%, g G W^, where V£, 
and Wfr are spaces of vector-valued functions, while W 2 
and Wfo are spaces of scalar-valued functions. The form 
takes four arguments, where three of the arguments are 
coefficients, 

a : V,] x Wl x W 2 x Wl -> R. (15) 
Hence, r = 1 and n = 3. 

Example 2.4. The H 1 ^) norm of the error e — u — Uh 
squared is 

a{;u,u h ) = I (u - u h ) 2 + \V(u - u h )\ 2 dx. (16) 
Jn 

The form takes two arguments and both are coefficients, 

a:Wlx W 2 R. (17) 
Hence, r = and n = 2. 

Defining variational forms for coupled PDEs can be per- 
formed in two ways in the above described framework. One 
approach is to couple the variational forms on the linear 
algebra level, using block vectors and block matrices and 
defining one form for each block. Alternatively, a single 
form for the coupled system may be defined using mixed 
finite elements. 



3 Finite Element Assembly 



To see how to compute the tensor A by summing the local 



contributions from each cell K, we let w? K 



\V K \ denote 



The standard algorithm for computing the global sparse the dimension of the local finite element space on K for 



tensor A is known a s assembly, see iZienkiewicz et al 
( 1967t ): iHughesI (|l987tl . By this algorithm, the tensor A 
may be computed by assembling (summing) the contribu- 
tions from the local entities of a finite element mesh. To 
express this algorithm for assembly of the global sparse 
tensor A for a general finite element variational form of 
rank r, we introduce the following notation and assump- 
tions. 

Let T = {K} be a set of disjoint cells (a triangulation 
or tesselation) partitioning the domain D, — UkctK- Fur- 
ther, let d e T denote the set of exterior facets (the set of 
cell facets on the boundary dil), and let d{T denote the set 
of interior facets (the set of cell facets not on the boundary 
dfl). For each discrete function space V£, j = 1, 2, . . . , r, 
we assume that the global basis is obtained by 

patching together local function spaces V K on each cell K 
as determined by a local-to-global mapping. 

We shall further assume that the variational form <j6j) 
may be expressed as a sum of integrals over the cells T, the 
exterior facets d e T and the interior facets d{T '. We shall 
allow integrals expressed on disjoint subsets T — ^l r =l Tk, 
d e T = ^l" =1 deTk and d{T = \J$L x di% respectively. 

We thus assume that the form a is given by 

a(v 1 ,...,v r ;w 1 ,...,w n ) = 

EE/ I k{ v lT--, v r]Wii...W n )&X 



E E / Ik(vi,---,Vr;W!,. 

k=l S<Ed a T k S 



,w n ) ds 



(18) 



E 



E 

k=l SediTk 



IKvx, . . . ,v r ;wi, . . . ,w n ) ds. 



We refer to an integral over a cell if as a cell integral, an 
integral over an exterior facet S as an exterior facet in- 
tegral (typically used to implement Neumann and Robin 
type boundary conditions), and to an integral over an in- 
terior facet S as an interior facet integral (typically used 
in discontinuous Galerkin methods). 

For simplicity, we consider here initially assembly of the 
global sparse tensor A corresponding to a form a given by 
a single integral over all cells T, and later extend to the 
general case where we must also account for contributions 
from several cell integrals, interior facet integrals and ex- 
terior facet integrals. 

We thus consider the form 



a(vi,...,v r ;wi,...,w n ) 



E 



I c {vi, . . . ,v r ;wi, . . . ,w n )dx, 



K 



for which the global sparse tensor A is given by 



E 



K 



.,4>l :wi,...,w n )dx. 



(19) 



(20) 



the jth primary argument Vj £ V fl for j = 1,2,.. 
Furthermore, let 



K 



[l,n j 



[l,N j ] 



(21) 



denote the local-to-global mapping for V£, that is, on any 
given K e T, the mapping l j k maps the number of a local 
degree of freedom (or, equivalently, local basis function) to 
the number of the corresponding global degree of freedom 
(or, equivalently, global basis function). We then define 
for each K 6 7" the collective local-to-global mapping lk '■ 
1 K -> I by 



L K{i) = {i K (il), i K {i2), ■ ■ ■ , tR^r)) Vi Elj 



K ■ 



(22) 



where Xk is the index set 



= (1,1, •••,2),..., (n K ,n 2 K ,..., rf K )}. 

(23) 

Furthermore, for each Vl we let {</>f r '"'}™ = fi denote the 
restriction to an element K of the subset of the basis 
{(t>{}^li C V K of V 3 h supported on K. 

We may now compute A by summing the contributions 
from the local cells, 

Ai= J2 I i c (ri 1 ,---,tf r ;w 1 ,...,w n )dx 

= E / / ^g ) - i (ii)'---'^ r )- i w ;u;i '---' u '" )dx 

KeTi JK 

= E K- K \w 

(24) 



where A K is the local cell tensor on cell K (the "element 
stiffness matrix" ) , given by 



A 



K 



r 



K 



iK,l 



, 4>f ,r ;w 1 ,...,w n )dx, 



(25) 



and where % denotes the set of cells on which all basis 
functions <p\ , <fq , . . . , 4>i r are supported. Similarly, we may 
sum the local contributions from the exterior and interior 
facets in the form of local exterior facet tensors and inte- 
rior facet tensors. 

In Algorithm Q] we present a general algorithm for as- 
sembling the contributions from the local cell, exterior 
facet and interior facet tensors into a global sparse ten- 
sor. In all cases, we iterate over all entities (cells, exterior 
or interior facets), compute the local cell tensor A K (or 
exterior/interior facet tensor A s ) and add it to the global 
sparse tensor as determined by the local-to-global map- 
ping, see Figure [TJ 



4(1) 4(2) 4(3) 



1 2 3 




Figure 1: Adding the entries of a cell tensor A K to the 
global tensor A using the local-to-global mapping lk, il- 
lustrated here for a rank two tensor (a matrix). 



Algorithm 1 Assembling the global tensor A from the lo- 
cal contributions on all cells, exterior and interior facets. 
For assembly over exterior facets, K(S) refers to the cell 
K G T incident with the exterior facet S, and for assembly 
over interior facets, K (S) refers to the "macro cell" con- 
sisting of the pair of cells K + and K~ incident with the 

interior facet S. 

A = 

(i) Assemble contributions from all cells 
for each K G T 



for j = 1,2, 



Tabulate the local-to-global mapping 



K 



for j = 1,2, 



Extract the values of mj, on K 



4 Software Framework for Finite Element Assembly 



In a finite element application code, typical input from the 
user is the variational (weak) form of a PDE, a choice of 
finite elements, a geometry represented by a mesh, and 
user-defined functions that appear as coefficients in the 
variational form. For a linear PDE, the typical solution 
procedure consists of first assembling a (sparse) linear sys- 
tem AU = b from given user input and then solving that 
linear system to obtain the degrees of freedom U for the 
discrete finite clement approximation Uh of the exact so- 
lution u of the PDE. Even when the solution procedure 
is more involved, as for a nonlinear problem requiring an 
iterative procedure, each iteration may involve assembling 
matrices and vectors. It is therefore clear that the assem- 
bly of matrices and vectors (or in general tensors) is an 
important task for any finite element software framework. 
We refer to the software component responsible for assem- 
bling a global tensor from given user input consisting of a 
variational form, finite element function spaces, mesh and 
coefficients as the Assembler. 

As demonstrated in Algorithm[TJ the Assembler needs to 
iterate over the cells in the mesh, tabulate degree of free- 
dom mappings, extract local values of coefficients, com- 
pute the local element tensor, and add each element ten- 
sor to the global tensor which is the final output. Thus, 
the Assembler is a software component where many other 
components are combined. It is therefore important that 
the software components on which the Assembler depends 
have well-defined interfaces. We discuss some issues re- 
lating to the design of these software components below 
and then demonstrate how these software components to- 
gether with the Assembler may be combined into a general 
software framework for finite clement assembly. 

4.1 Variational Forms 

Implementations of discrete variational forms in a general 
finite clement library usually consist of programming ex- 



Take < k < n c such that K G Tk 
Tabulate the cell tensor A K for If. 
Add Af to A L i K{ii y L 2 K{l2) ^ K{ir) for i G I K 

(ii) Assemble contributions from all exterior facets 
for each S G d e T 



for j = l,2,...,r: 

Tabulate the local-to-global mapping v 

for j = 1,2,.,., n: 

Extract the values of Wj on K(S) 

Take < k < n e such that S G d e Tk 
Tabulate the exterior facet tensor A s for It 



K{S) 



Add Af to A,i 



K(S)( l l)' t K(S)( i 2). 



K(S 



(i r ) for i G I K (s) 



(iii) Assemble contributions from all interior facets 
for each S G d{T 



for j = 1,2,..., r: 

Tabulate the local-to-global mapping v 

for j = l,2,...,n: 

Extract the values of Wj on K(S) 



K(S) 



Take < k < n.i such that S G d{Tk 



Tabulate the interior facet tensor A s for II 



Add Af to A,i 



K(S) (* 1 )' t K(S) (* 2 )> 



K(S) 



(ir) 



for i G Ik 



(S) 



pressions for the integrands If, If, If (see Equation (TT5)) ). 
eventually writing a quadrature loop and a loop over el- 
ement matrix indice s depending on t he ab st raction level 
of the library (see Bangerth et al. ( 2007t h iLangtangen 
( 2003tJ )). An alternative approach is to apply exact in- 
tegration instead of quadrature. In either case, the result 
of this computation may be communicated through the 



UFC interface. 

The motivation behind the UFC interface is to separate 
the implementation of the form from other details of the 
assembly such as the mesh and the linear algebra libraries 
in use. 

In the FEniCS finite element software framework, a 
high-level form language embedded in Python is used to 
define the variational form and finite elements. This re- 
duces the distance from the mathematical formulation of 
a PDE to an implementation of a PDE solver, removes 
tedious and error-prone tasks such as coding PDE-specific 
assembly loops, and enables rapid prototyping of new mod- 
els and methods. To retain computational efficiency, we 
generate efficient low-level code from the abstract form 
description, using exact integration where possible. Code 
generation adds another complexity layer to the software, 
and it becomes even more important to keep a clear separa- 
tion between software components such that the interface 
between generated code and library code is well defined. 
This is achieved by generating implementations of the UFC 
interface. 

4.2 Mesh Libraries 

Many different representations of computational meshes 
exist. Typically, each finite element library provides its 
own internal implementation of a computational mesh. We 
do not wish to tie the UFC interface to one particular mesh 
representation or one particular library. Still, several op- 
erations like the element tensor computation depends on 
local mesh data. For this reason, the UFC interface pro- 
vides a low-level data structure to communicate single cell 
data. In addition, a small data structure is used to com- 
municate global mesh dimensions which are necessary for 
computing the local-to-global mapping. Assemblers im- 
plemented on top of the UFC interface must therefore be 
able to copy/translate cell data from the mesh library be- 
ing used to the UFC data structure (involving a mini- 
mal overhead). This makes it possible to achieve sepa- 
ration between the mesh representation and the element 
tensor computation. The Assembler implemented in FEn- 
iCS (as part of D OLFIN) is im plemented for one particular 
mesh format, see lLogd (120081) . but an Assembler compo- 



nent could easily be written for other mesh libra ries like 
the PETSc Sieve, see lKarpeev and Kneplevl (|2005l) and the 
DUNE-Grid interface ((iBastian et all <|2007allbh ). 



4.3 Linear Algebra Libraries 

It is desirable to reuse existing high-performance linear 
algebra libraries like PETSc and Trilinos. It is therefore 
important that the Assembler is able to assemble element 
tensors into global matrices and vectors implemented by 
external libraries. Aggregation into matrix and vector data 
structures are fairly similar operations, typically consisting 
of passing array pointers to existing functions in the linear 
algebra libraries. By implementing a common interface for 
assembly into tensors of arbitrary rank, the same assembly 



routine can be reused for any linear algebra library without 
changes. This avoids duplication of assembly code, and 
one may easily change the output format of the assembly 
procedure. The details of these interfaces are beyond the 
scope of the current paper. At the time of writing, we 
have written assembly routines with support for matrices 
and vectors from Epetra (Trilinos), PETSc, PyCC, and 
uBLAS in addition to scalars. 

With components available for finite element variational 
forms, mesh representation, and linear algebra, we may 
use the UFC interface to combine these components to 
build a PSE for partial differential equations. The central 
component of this PSE is the Assembler. As illustrated in 
Figure [2j the Assembler takes as input a variational form, 
communicated through the UFC interface, a mesh and a 
set of functions (the coefficients), and assembles a tensor. 
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Figure 2: Assembling a tensor from a given UFC, mesh 
and functions (coefficients). 



4.4 High-level Interfaces 

In FEniCS, we have additional application-level abstrac- 
tions for expressing variational forms, meshes, functions 
and linear algebra objects to achieve a consistent high-level 
user-interface. The generation of the UFC code may then 
be hidden from the user, who just provides a high-level 
description of the form. The PSE may then automatically 
generate the UFC at run-time, functioning as just-in-time 
(JIT) compiler, and call the Assembler with the generated 
UFC. Below, we demonstrate how this may be done in the 
Python interface of DOLFIN. The user here defines a finite 
element function space, and a pair of bilinear and linear 
forms a(v, u) = J n Vv • Vu + vu dx and L = J n vf dx, from 
which a matrix and vector may be assembled by calls to 
the function assemble. A linear system solver may then 
be invoked to compute the degrees of freedom U of the 
solution. 



element = FiniteElementC'CG" , "triangle" 

v = TestFunctionCelement) 

u = TrialFunction(element) 

f = FunctionCelement , mesh, 100.0) 

a = dot(gradCv), gradCu))*dx + v*u*dx 
L = v*f*dx 

A = assemble Ca, mesh) 
b = assemble CL, mesh) 



U = solve (A, b) 



5 The UFC Interface 



While FEniCS provides an integrated environment, in- 
cluding the PSE DOLFIN and the two form compilers FFC 
(FEniCS Form Compiler) and SFC (SyFi Form Compiler), 
the fact that these components comply with the UFC inter- 
face means that they may also be used interchangeably in 
heterogeneous environments together with other libraries 
(that implement or use the UFC interface). This is il- 
lustrated in the flow diagram of Figure [3] where alternate 
routes from mathematical description to matrix assembly 
are demonstrated. 
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Figure 5: UML diagram of the UFC class relations 
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Figure 3: Alternate routes from mathematical description 
to matrix assembly enabled by the UFC interface (Note 
that the Diffpack example is fictional). 



In Figure [31 we have also included another interface UFL 
(Unified Form Language) which provides a unified way to 
express finite element variational forms. The UFL inter- 
face is currently in development. Together, UFL and UFC 
provide a unified interface for the input and output of form 
compilers, see Figure H]. 
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Figure 4: An abstract definition (UFL) of a finite element 
variational form is given as input to a form compiler, which 
generates UFC code as output. 



The UFC interface consists of a small collection of ab- 
stract C++ classes that represent common components 
for assembling tensors using the finite element method. 
These classes are accompanied by a well documented 
(jAlnass et all (|2008tl ) set of conventions for numbering of 
cell data and other arrays. We have strived to make the 
classes as simple as possible while not sacrificing gener- 
ality or efficiency. Data is passed as plain C arrays for 
efficiency and minimal dependencies. Most functions are 
declared const, reflecting that the operations they repre- 
sent should not change the outcome of future operations!! 
Other initialization of implementation-specific data should 
ideally be performed in constructors. 

One can ask why the UFC interface consists of classes 
and not plain functions. There are three reasons for this. 
First, we want to handle each form as a self contained 
"black box" , which can be passed around easily in an ap- 
plication. Many functions belong together conceptually, 
thus making it natural to collect them in a class "names- 
pace" . Second, we need multiple versions of each function 
in the software representation of variational forms, in par- 
ticular to represent multiple variational forms and multi- 
ple finite element function spaces. This is best achieved by 
making each such function a member function of a class 
and having multiple instances of that class. Third, UFC 
function implementations may need access to stored data, 
and with a plain function-based interface these data would 
then need to be global variables. In particular, when exist- 
ing libraries or applications want to implement the UFC 
interface, it may be necessary for the subclasses of UFC 
classes to inherit from existing classes or to have pointers 
to other objects. 

5.1 Class Relations 

Figure ([5]) shows all the classes and their relations. The 
classes mesh, cell, and function provide the means for 



The exceptions are the functions to initialize a dofjnap. 



communicating mesh and coefficient function data as ar- 
guments. Each argument of the form (both primary argu- 
ments and coefficients) is represented by a f inite_element 
and dof _map object. The integrals are represented by one 
of the classes cell_integral, exterior jfacet_integral, 
or interior j£acet_integral. An object of the class form 
gives access to all other objects in a particular implemen- 
tation. In this paper, we will not describe all the functions 
of these classes in de tail. A complete sp ecification can be 
found in the manual ( Alnaes et al. ( 2008t )V 

At the core of UFC is the class form, which 
represents the general variational form a of Equa- 
tion (|18p . Subclasses of form must implement 
factory functions which may be called to cre- 
ate cell_integral, exterior jfacet_integral and 
interiorJ: acet_integral objects. These objects in turn 
know how to compute their respective contribution from 
a cell or facet during assembly. A code fragment from the 
form class declaration is shown below. 



class form 
{ 

public : 



/// Create cell integral on sub domain i 
virtual cell_integral* 
create_cell_integral (unsigned int i) 
const = 0; 

/// Create exterior facet integral on sub domain i 
virtual exterior_f acet_integral* 
create_exterior_facet_integral (unsigned int i) 
const = 0; 

/// Create interior facet integral on sub domain i 
virtual interior_f acet_integral* 
create_interior_facet_integral (unsigned int i) 
const = 0; 



The form class also specifies functions for creating 
f inite_element and dof _map objects for the finite clement 
function spaces {V£}j =1 , {W^}™ =1 of the variational form. 
The f inite_element object provides functionality such as 
evaluation of degrees of freedom and evaluation of basis 
functions and their derivatives. The dof_map object pro- 
vides functionality such as tabulating the local-to-global 
mapping of degrees of freedom on a single element, as well 
as tabulation of subsets associated with particular mesh 
entities, used to apply Dirichlet boundary conditions and 
build connectivity information. 

Both the f inite_element and dof_map classes can rep- 
resent mixed elements, in which case it is possible to ob- 
tain f inite_element and dof_map objects for each sub- 
element in a hierarchical manner. Vector elements com- 
posed of scalar elements are in this context seen as special 
cases of mixed elements where all sub-elements are equal. 
Thus, e.g., from a dof _map representing a P2 — Pi Taylor- 
Hood element, it is possible to extract one dof _map for the 
quadratic vector element and one dof _map for the linear 
scalar element. From the vector element, a dof_map for 
the quadratic scalar element of each vector component can 



be obtained. This can be used to access subcomponents 
from the solution of a mixed system. 

5.2 Stages in the Assembly Algorithm 



enum shape {interval, triangle, quadrilateral, 
tetrahedron, hexahedron}; 

class cell { 
public : 

shape cell_shape; 

unsigned int topological_dimension; 
unsigned int geometric_dimension; 

/// Array of global indices for the mesh entities of the cell 
unsigned int** entity_indices ; 

/// Array of coordinates for the vertices of the cell 
double** coordinates; 



Figure 6: 
data. 



Data structure for communicating single cell 



Next, we focus on a few key parts of the interface and 
explain how these can be used to implement the assembly 
algorithm (Algorithm [IJ . This algorithm consists of three 
stages: (i) assembling the contributions from all cells, (ii) 
assembling the contributions from all exterior facets, and 
(iii) assembling the contributions from all interior facets. 

Each of the three assembly stages (i)-(iii) of Algorithm!]] 
is further composed of five steps. In the first step, the 
polygon K is fetched from the mesh, typically implemented 
by filling a cell structure (see Figure |6]) with coordinate 
data and global numbering of the mesh entities in the cell. 
This step depends on the specific mesh being used. 

Secondly, the local-to-global mapping of degrees of free- 
dom is tabulated for each of the function spaces. That 
is, for each of the discrete finite element spaces on K, we 
tabulate (or possibly compute) the global indices for the 
degrees of freedom on {V^}£ =1 and 

The class dof _map represents the mapping between local 
and global degrees of freedom for a finite element space. 
A dof_map is initialized with global mesh dimensions by 
calling the function init_mesh (const meshfe m) . If this 
function returns true, the dof _map should be additionally 
initialized by calling the function init_cell (const meshfe 
m, const cellfe c) for each cell in the global mesh, fol- 
lowed by init_cellJ:inalize after the last cell. After 
the initialization stage, the mapping may be tabulated at 
a given cell by calling a function with the following signa- 
ture. 



void dof _map :: tabulate_dofs (unsigned int* dofs, 
const mesh& m, 
const cell& c) const 



Here, unsigned int* dofs is a pointer to the first el- 
ement of an array of unsigned integers that will be filled 
with the local-to-global mapping on the current cell during 
the function call. 



In the third step of each stage of Algorithm [TJ we may 
use the tabulated local-to-global mapping to interpolate 
(extract) the local values of any of the coefficients in 

If a coefficient Wj is not given as a linear combination 
of basis functions for , it must at this step be interpo- 
lated into Wl , using the interpolant defined by the degrees 
of freedom of Wl (for example point evaluation at a set 
of nodal points). In this case, the coefficient function is 
passed as an implementation of the function interface (a 
simple functor) to the function evaluate_dof s. 

/// Evaluate linear functionals for all dofs on the function f 

void f inite_element : :evaluate_dof s(double *values, 

const f unction& f , 
const cell& c) const 



In the fourth step, the local element tensor contribu- 
tions (cell or exterior/interior facet tensors) are computed. 
This is done by a call to the function tabulate_tensor, 
illustrated below for a cell integral. 

void cell_integral :: tabulate_tensor (double* A, 

const double * const * w, 
const cell& c) const 



Similarly, one may evaluate interior and exterior facet con- 
tributions using slightly different function signatures. 

Finally, at the fifth step, the local element tensor contri- 
butions are added to the global tensor, using the local- 
to-global mappings previously obtained by calls to the 
tabulate_dof s function. This is a simple operation that 
depends on the linear algebra library in use. 



6 Examples 



In this section, we demonstrate how UFC is used in prac- 
tice in DOLFIN, FFC, and SFC. First, we show a part 
of the assembly algorithm (Algorithm [TJ) as implemented 
in DOLFIN. We then show examples of input to the form 
compilers FFC and SFC as well as part of the correspond- 
ing UFC code generated as output. Examples include Pois- 
son's equation and linear convection (see Example I2.2|) . 

6.1 An Example UFC Assembler 

To demonstrate how one may implement an assembler 
based on the UFC interface, we provide here a (somewhat 
simplified) excerpt from the DOLFIN assembler!! 

for (Celllterator cell (mesh); Icell.endO; ++cell) 
{ 

uf c .update (* cell) ; 

for (uint i = 0; i < uf c . f orm . rankO ; i++) 

dof _map_set [i] . tabulate_dof s (uf c . dof s [i] , *cell) ; 

for (uint i = 0; i < coef f icients . size () ; i++) 

3 The ufc object is here an instance of a simple DOLFIN class 
that holds pointers to arrays and UFC container classes, such as the 
array A and cell data ufc: :cell, needed to communicate through the 
UFC interface. 



coef f icients [i] -interpolate (ufc . w [i] , ufc. cell, 

*uf c . coef f icient_elements [i] , 
♦cell) ; 

integral->tabulate_tensor (uf c . A, ufc.w, ufc. cell); 
A.addCufc.A, uf c . local_dimensions , ufc. dofs); 

} 



The outer loop iterates over all cells in a given mesh. For 
each cell, a ufc: : cell is updated and the local-to-global 
mapping is constructed. We then interpolate all the form 
coefficients on the cell and compute the element tensor. 
At the end of the iteration, the local-to-global mapping is 
used to add the local tensor to the global tensor. 

6.2 FFC Examples 

The form compiler FFC provides a simple language for 
specification of variational forms, which may be entered ei- 
ther directly in Python or in text files given to the compiler 
on the command-line. For each variational form given as 
input, FFC generates UFC-compliant C++ code for eval- 
uation of the corresponding element tensor (s). 



Poisson's Equation 

We consider the following input file to FFC for Poisson's 
equation. 



element = FiniteElement("Lagraiige" , 


"triangle" , 1) 


v = TestFunction(element) 




u = TrialFunction(element) 




f = Function(element) 




a = dot(gradCv), grad(u))*dx 




L = v*f*dx 





Here, two forms a (bilinear) and L (linear) are defined. 
Both the test and trial spaces are spanned by linear La- 
grange elements on triangles in two dimensions. When 
compiling this code using FFC, a C++ header file is cre- 
ated, containing UFC code that may be used to assemble 
the global sparse stiffness matrix and load vector. Below, 
we present the code generated for evaluation of the element 
stiffness matrix for the bilinear form a. 



virtual void tabulate_tensor (double* A, 

const double * const * w, 
const ufc::cell& c) const 

{ 

// Extract vertex coordinates 

const double * const * x = c . coordinates ; 

// Compute Jacobian of affine map from reference cell 
const double J_00 = x[l][0] - x[0] [0] ; 
const double J_01 = x [2] [0] - x[0] [0] ; 
const double J_10 = x[l][l] - x[0] [1] ; 
const double J_ll = x[2][l] - x[0] [1] ; 

// Compute determinant of Jacobian 
double detJ = J_00*J_11 - J_01*J_10; 

// Compute inverse of Jacobian 
const double Jinv_00 = J_ll / detJ; 
const double Jinv_01 = -J_01 / det J; 
const double Jinv_10 = -J_10 / detJ; 
const double Jinv_ll = J_00 / detJ; 

// Set scale factor 

const double det = std : : abs (det J) ; 



// Compute geometry tensors 

const double G0_0_0 = det* (Jinv_00*Jinv_00 + Jinv_01*Jinv_01) 
const double G0_0_1 = det* (Jinv_00*Jinv_10 + Jinv_01*Jinv_ll) 
const double G0_1_0 = det* (Jinv_10*Jinv_00 + Jinv_ll*Jinv_01) 
const double G0_1_1 = det* ( Jinv_10*Jinv_10 + Jinv_ll*Jinv_ll) 



// Compute element tensor 

A[0] = 0.5*G0_0_0 + 

A[l] = -0.5*G0_0_0 

A [2] = -0.5*G0_0_1 

A [3] = -0.5*GCLCL0 

A [4] = 0.5*G0_0_0; 

A [5] = 0.5*G0_0_lj 

A [6] = -0.5*G0_1_0 

A [7] = 0.5*G0_1_0; 

A[8] = 0.5*G0_1_1; 



0.5*G0_0_1 + 0.5*G0_1. 

■ 0.5*G0_1_0; 
- 0.5*G0_1_1; 

■ 0.5*G0_0_1; 



0.5*G0_1_1; 



+ 0.5*G0_1_1; 



In FFC, an element tensor contribution is computed as 
a tensor contraction between a geometry tensor varying 
from cell to cell, and a geometry indep e ndent tenso r on a 
reference element, see iKirbv and Loggl (|2006l 120071 ). For 
simple forms, like the one under discussion, the main work 
is then to construct the geometry tensor, related to the 
geometrical mapping between the reference element and 
physical element. 

Having computed the element tensor, one needs to com- 
pute the local-to-global mapping in order to know where 
to insert the local contributions in the global tensor. This 
mapping may be obtained by calling the member func- 
tion tabulate_dof s of the class uf c: :dof_map. FFC uses 
an implicit ordering scheme, based on the indices of the 
topological entities in the mesh. This information may be 
extracted from the ufc: :cell attribute entity_indices. 



virtual void tabulate_dofs (unsigned int* dofs, 
const ufc::mesh& m, 
const ufc::cell& c) const 



< 



dofs[0] = 
dofs[l] = 
dofs [2] = 



. entity_indices [0] [0] 
. entity_indices [0] [1] 
. entity_indices [0] [2] 



For Lagrange elements on triangles, each degree of free- 
dom is associated with a global vertex. Hence, FFC con- 
structs the mapping by picking the corresponding global 
vertex number for each degree of freedom. 



Linear Convection 

Consider the variational form in Example 
file to FFC reads as follows. 



The input 



vector_element = VectorElement("Lagrange" 


"triangle" 


1) 


scalar_element = FiniteElementC'Lagrange" 


"triangle" 


1) 


v = TestFunction(vector_element) 






u = TrialFunction(vector_element) 






w = Function(vector_element) 






rho = Function(scalar_element) 






a = rho*v [i] *w [j] *u [i] . dx ( j ) *dx 







The code generated for the tabulate_tensor function is 
presented below. Computations involving coefficients are 
performed by interpolating the functions w and p on the 
cell under consideration. These values are stored in the 



array w below. For clarity, some code has been omitted in 
this example. 



virtual void tabulate_tensor (double* A, 

const double * const * w, 
const ufc::cell& c) const 

{ 

// Extract vertex coordinates and compute Jacobian etc 
// as in previous example 



// Compute coefficients 

const double cl_0_0_0 = w[l] [0] 

const double cl_0_0_l = u[l] [1] 

const double cl_0_0_2 = w[l] [2] 

const double c0_0_l_0 = w[0] [0] 

const double c0_0_l_5 = w[0] [5]; 
// Compute geometry tensors 

const double G0_0_0_0_0 = det*cl. 

const double G0_0_0_1_0 = det*cl. 

const double G0_0_1_0_0 = det*cl. 



0_0_0*c0_0_l_0*Jinv_00 
0_0_0*cO_0_l_0*Jinv_10 
0_0_0*c0_0_l_l*Jinv_00 



// Compute element tensor 
A[0] = -0.05*G0_0_0_0_0 - 
A[l] = 0.05*G0_CLCLCL0 + 
A [2] = 0.05*G0_0_0_1_0 + 



The local-to-global mapping for the space of piecewise 
linear vectors is computed by associating two values with 
each vertex. The code generated for tabulate_dof s is 
presented below. 



virtual void tabulate_dof s (unsigned int* dofs, 
const ufc: :mesh& m, 
const ufc::cell& c) const 



{ 



dofs[0] = c . entity_indices [0] [0] 

dofs[l] = c . entity_indices [0] [1] 

dofs [2] = c.entity_indices[0] [2] 
unsigned int offset = m.num_entities [0] ; 

dofs [3] = offset + c . entity_indices [0] [0] 

dofs[4] = offset + c . entity_indices [0] [1] 

dofs [5] = offset + c . entity_indices [0] [2] 



FFC generates code for arbitrary multilinear forms and 
currently supports arbitrary degree continuous Lagrange 
elements, discontinuous elements, RT elements, BDM ele- 
ments, BDFM elements and Nedelec elements in two and 
three space dimensions. 

6.3 SFC Examples 

SFC is another form compiler producing UFC code, in 
which the user defines variational forms in Python us ing a 
symbolic engine based on GiNaC ([Bauer et al. | d2006h ). It 
has a slightly different feature set than FFC, such as us- 
ing symbolic differentiation to automatically compute the 
Jacobi matrix of a nonlinear form. The resulting low-level 
UFC code is very similar. 

Power-law Viscosity 

Example 12.31 is specified in SFC as follows. 



vector_element = VectorElementC'Lagrange" , "triangle", 1) 
scalar_element = FiniteElementC'DG" , "triangle", 0) 

v = TestFunction(vector_element) 
w = Function(vector_element) 



mu = FunctionCscalar_element) 
rho = Function(scalar_element) 

def power_law(v, w, mu, rho, itg) : 
q = 0.3 

GinvT = itg.GinvTO 
Dw = grad(w, GinvT) 
Dv = gradCv, GinvT) 
wDw = dot(w, Dw) 

return rho*dot (wDw, v) + mu*inner(Dw,Dw)**q * inner(Dw,Dv) 

F_form = FormCbasisf unctions = [v] , 

coefficients = [w, mu, rho]) 
F_f orm . add_cell_ integral (power_law) 
J_form = Jacobi (F_f orm) 



The syntax for denning elements and arguments is the 
same as in FFC, but the integrand is specified in a slightly 
different syntasQ This code also computes the form cor- 
responding to the Jacobian matrix, using symbolic differ- 
entiation. The generated code for computing the Jacobian 
matrix is in this case more complicated but it implements 
the UFC interface in the same manner as in the previous 
examples. 

In the current implementation, SFC explicitly constructs 
a local-to-global mapping at run-time. In this case, with 
Lagrange elements, the global coordinates identify the de- 
grees of freedom. The UFC interface supports constructing 
the local-to-global mapping through the initjnesh and 
init_cell methods of ufc: : dof jnap. Below, we present 
the code generated for init_cell (where we use additional 
structures of type Ptv (point) for representing degrees of 
freedom and the container Dof_Ptv dof for building the 
local-to-global mapping) . 

void dof _map_2D :: init_cell (const ufc::mesh& m, const ufc::cell& c) 
{ 

// coordinates 

double xO = c . coordinates [0] [0] ; double yO = c . coordinates [0] [1] ; 

double xl = c . coordinates [1] [0] ; double yl = c . coordinates [1] [1] ; 

double x2 = c . coordinates [2] [0] ; double y2 = c . coordinates [2] [1] ; 

// affine map 

double GOO = xl - xO; 

double G01 = x2 - xO; 

double G10 = yl - yO; 
double Gil = y2 - yO; 

unsigned int element = c . entity_indices [2] [0] ; 

double dof0[2] = { xO, yO }; 

Ptv pdof0(2, dofO) ; 

dof . insert_dof (element , 0, pdof 0) ; 

double dofl[2] = { GOO+xO, yO+GlO }; 

Ptv pdof 1(2, dof 1) ; 

dof . insert_dof (element , 1, pdof 1) ; 

double dof 2 [2] = { xO+GOl, Gll+yO }; 

Ptv pdof 2 (2, dof 2) ; 

dof . insert_dof (element , 2, pdof 2) ; 

} 



The dofjnap class is only responsible for the unique- 
ness of the local-to-global mapping. Possible renumbering 
strategies may be imposed by the assembler, for example 
to minimize communication when assembling in parallel. 

4 The variable itg is an integral object containing information 
about the mapping between physical coordinates and the reference 
element. 



7 Discussion 



We have used (generated) UFC for many applications, in- 
cluding Poisson's equation; convection-diffusion-reaction 
equations; continuum equations for linear elasticity, hy- 
perelasticity, and plasticity; the incompressible Navier- 
Stokes equations; and mixed formulations for the Hodge 
Laplacian. The types of finite elements involved in- 
clude standard continuous Lagrange elements of arbi- 
trary order, discontinuous Galcrkin formulations, BDM el- 
ements, Raviart-Thomas elements, Crouzeix-Raviart ele- 
ments, and Nedelec elements. 

The form compilers FFC and SFC are UFC compli- 
ant, both generating efficient UFC code from an abstract 
problem definition. Assemblers have been implemented in 
DOLFIN and PyCC, using the DOLFIN mesh represen- 
tation, and together covering linear algebra formats from 
PETSc, Trilinos (Epetra), uBLAS, and PyCC. Parallel as- 
sembly is supported in the current development version 
of DOLFIN, without requiring any modifications to UFC 
since it operates on an element level. Altogether, this 
demonstrates that the UFC interface is flexible both in 
terms of the applications and finite element formulations 
it covers, and in terms of its interoperability with existing 
libraries. 

One of the main limitations in the current version of the 
UFC interface (vl.l) is the assumption of a homogeneous 
mesh, that is, only one cell shape is allowed throughout 
the mesh. Thus, although mesh ordering conventions have 
been defined for the interval, triangle, tetrahedron, quadri- 
lateral, and hexahedron, only one type of shape can be 
used at any time. Also, higher order (non-affinely mapped) 
meshes are not supported in the current version of the in- 
terface. Another limitation is that only one fixed finite el- 
ement space can be chosen for each argument of the form, 
which excludes p-refinement (increasing the element or- 
der in a subset of the cells). All these limitations may 
be removed in future versions of UFC, and we encourage 
interested developers to make contact to address these lim- 
itations. 

UFC provides a unified interface for code generated as 
output by form compilers such as FFC and SFC. Similarly, 
we are currently working on a specification for a unified 
form language (UFL) to function as a common input to 
form compilers. Currently, both FFC and SFC provide 
(different) form languages for easy specification of varia- 
tional forms in a high-level syntax. With a unified form 
language, a user may specify a variational form in that 
language and assemble the corresponding discrete operator 
(tensor), independently of the components being used to 
generate the UFC code from the UFL, and independently 
of the components being used to assemble the tensor from 
the UFC form. 



8 Conclusion 



We have presented a general framework for assembly of 
finite element variational forms. Based on this framework, 
we have then extracted an interface (UFC) that may be 
used to provide a communication layer between general- 
purpose and problem-specific code for assembly of finite 
clement variational forms. 

The interface makes minimal assumptions on the type of 
problem being solved and the data structures involved. For 
example, the discrete variational form may in general be 
multilinear and hence assemble into a tensor of arbitrary 
rank. The basic data structures used to pass data through 
the interface are composed of plain C arrays. The minimal 
set of assumptions on problem and data structures enables 
application of the interface to a wide range of variational 
forms and a large collection of finite element libraries. 

We have used the UFC interface in the implementation 
of the FEniCS suite of finite element tools. In a simple 
Python script, one may define a variational form and a 
mesh, and assemble the corresponding global sparse ma- 
trix (or vector). When doing so, UFC code is generated 
by either of the form compilers FFC or SFC, and passed 
to the UFC-compatible assembler of the general-purpose 
finite element library DOLFIN. 

We encourage developers of finite element software to 
use the UFC interface in their libraries. By doing so, those 
libraries may directly take advantage of the form compil- 
ers FFC or SFC to specify finite clement problems. More- 
over, one can think of already existing specifications of 
complicated finite element problems that via UFC can be 
combined with other libraries than the specifications were 
originally written for. We have tried to make minimal as- 
sumptions to make this possible. 

We believe that UFC itself and the ideas behind it con- 
stitute an important step towards greater flexibility in fi- 
nite element software. By code generation via tools like 
FFC and SFC, this flexibility may be retained also in com- 
bination with very high performance. 



REFERENCES 



Cactus. http://www. cactuscode. org/. 

COMSOL Multiphysics. http://www.comsol.com. 

FEniCS. http://www.fenics.org. 

Getfem++. http://home.gna.org/getfem/. 

Hypre. http://acts.nersc.gov/hypre/. 

Kaskade. 

http://www.zib.de/Numerik/numsoft/kaskade/. 

Trilinos. http:/ /software. sandia.gov/trilinos. 

Alnass, M., H.-P. Langtangen, A.Logg, K.-A. Mardal, and 
O. Skavhaug (2008). UFC Specification and User Man- 
ual, http //www. fenics.org/ufc/. 



Alnaes, M. S. and K.-A. Mardal (2008). SyFi: Symbolic 
finite elements, http://www.fenics.org/syfi/. 

Balay, S., K. Buschelman, W. D. Gropp, D. Kaushik, M. G. 
Knepley, L. C. Mclnnes, B. F. Smith, and H. Zhang 
(2008). PETSc. http://www.mcs.anl.gov/petsc/. 

Bangcrth, W., R. Hartmann, and G. Kanschat 
(2006). deal . II Differential Equations Analysis Library. 
http: //www. dealii. org/. 

Bangcrth, W., R. Hartmann, and G. Kanschat (2007). 
Deal. II — a general-purpose object-oriented finite ele- 
ment library. ACM Trans. Math. Softw. 33(4). 

Bastian, P., K. Birkcn, S. Lang, K. Johannsen, N. Neuss, 
H. Rcntz-Rcichcrt, and C. Wieners (1997). UG - a flexi- 
ble software toolbox for solving partial differential equa- 
tions. Computing and Visualization in Science 1, 27-40. 

Bastian, P., M. Blatt, A. Dedner, C. Engwer, R. Klofkorn, 
M. Ohlberger, and O. Sander (2007a). A generic grid 
interface for parallel and adaptive scientific computing. 
Part I: Abstract framework. Computing. Submitted. 
Available from the DUNE website. 

Bastian, P., M. Blatt, A. Dedner, C. Engwer, R. Klofkorn, 
M. Ohlberger, and O. Sander (2007b). A generic grid 
interface for parallel and adaptive scientific computing. 
Part II: Implementation and tests in DUNE. Computing . 
Submitted. Available from the DUNE website. 

Bauer, C, C. Dams, A. Frink, V. V. Kisil, R. Kreckel, 
A. Sheplyakov, and J. Vollinga (2006). GiNaC. 
http:/ /www. ginac.de. 

Ciarlct, P. G. (1978). The Finite Element Method for El- 
liptic Problems. North-Holland, Amsterdam, New York, 
Oxford. 

Hughes, T. J. R. (1987). The Finite Element Method: 
Linear Static and Dynamic Finite Element Analysis. 
Prentice-Hall. 

Karpeev, D. A. and M. G. Knepley (2005). Flexible repre- 
sentation of computational meshes, submitted to ACM 
Trans. Math. Softw.. 

Kirby, R. C. (2004). FIAT: A new paradigm for comput- 
ing finite element basis functions. ACM Trans. Math. 
Software 30, 502-516. 

Kirby, R. C. and A. Logg (2006). A compiler for varia- 
tional forms. ACM Transactions on Mathematical Soft- 
ware 32(3), 417-444. 

Kirby, R. C. and A. Logg (2007). Efficient compilation 
of a class of variational forms. A CM Transactions on 
Mathematical Software 33(3). 

Langtangen, H. P. (2003a). Computational Partial Dif- 
ferential Equations - Numerical Methods and Diffpack 



Programming (2nd ed.). Texts in Computational Sci- 
ence and Engineering, vol 1. Springer. 
http://www. diffpack. com. 

Langtangen, H. P. (2003b). Computational Partial Dif- 
ferential Equations - Numerical Methods and Diffpack 
Programming. Springer- Verlag. 2nd edition, 855 pages. 

Logg, A. (2007). Automating the finite element method. 
Arch. Comput. Methods Eng. 14(2), 93-138. 

Logg, A. (2008). Efficient representation of computational 
meshes. Submitted to International Journal of Compu- 
tational Science and Engineering. 

Logg, A. et al. (2008). FFC: FEniCS form compiler. 
http//www.fenics. org/ ffc/. 

Logg, A., G. Wells, J. Hoffman, J. Jansson, et al. 
DOLFIN: A general-purpose finite element library. 
http:/ / www.fenics. org/ dolfin/. 

Long, K. (2006). Sundance. 

http://software.sandia.gov/sundance/. 

Meerbergen, K. (2008). GLAS: Generic interface for Linear 
Algebra Software, http://glas.sourceforge.net/. 

Ring, J., H. P. Langtangen, and R. Bredesen (2008). 
Easyviz. Subpackage of SciTools: 
http://code.google.eom/p/scitools/. 

Skavhaug, O., M. S. Alnaes, K.-A. Mardal, G. Staff, 
A. Odegard, and et al (2008). PyCC. Software frame- 
work under development. 
http:/ / www. Simula, no /pyec. 

Thunc, M., E. Mossbcrg, P. Olsson, J. Rantakokko, 
K. Ahlander, and K. Otto (1997). Object-oriented con- 
struction of parallel PDE solvers. In E. Argc, A. M. 
Bruaset, and H. P. Langtangen (Eds.), Modern Software 
Tools for Scientific Computing, pp. 203-226. Birkhauser. 

Zienkiewicz, O. C, R. L. Taylor, and J. Z. Zhu (2005, hrst 
published in 1967). The Finite Element Method — Its 
Basis and Fundamentals, 6th edition. Elsevier. 



